Accurate Solution of Structured Least Squares Problems via Rank-Revealing Decompositions

Least squares problems minx ‖b−Ax‖2 where the matrix A ∈ Cm×n (m ≥ n) has some particular structure arise frequently in applications. Polynomial data fitting is a well-known instance of problems that yield highly structured matrices, but many other examples exist. Very often, structured matrices have huge condition numbers κ2(A) = ‖A‖2 ‖A†‖2 (A† is the Moore-Penrose pseudo-inverse ofA) and, therefore, standard algorithms fail to compute accurate minimum 2-norm solutions of least squares problems. In this work, we introduce a framework that allows us to compute minimum 2-norm solutions of many classes of structured least squares problems accurately, i.e., with errors ‖x̂0 − x0‖2/‖x0‖2 = O(u), where u is the unit roundoff, independently of the magnitude of κ2(A) for most vectors b. The cost of these accurate computations is O(n2m) flops, i.e., roughly the same cost as standard algorithms for least squares problems. The approach in this work relies in computing first an accurate rank-revealing decomposition of A, an idea that has been widely used in the last decades to compute, for structured ill-conditioned matrices, singular value decompositions, eigenvalues and eigenvectors in the Hermitian case, and solutions of linear systems with high relative accuracy. In order to prove that accurate solutions are computed, it is needed to develop a multiplicative perturbation theory of least squares problems. The results presented in this paper are valid in the case of both full rank and rank deficient problems and also in the case of underdetermined linear systems (m < n). Among other types of matrices, the new method applies to rectangular Cauchy, Vandermonde, and graded matrices and detailed numerical tests for Cauchy matrices are presented.